*!Version 1.0 29Aug13

*** 8/13 takes out additional weights, modified from sgnullsumbsamplez-1.0.ado
*** No longer constrains means to be the same across T AND C

*** sgnullsubsamplez.ado
*** Taken from version 3 of sgnullsubsample
*** adds 2 outcomes which keep zeros in for figures with zeros

** CHANGELOG ** for sgnullsubsample
*!Version 3.0.0 14Dec11
*!Version 2.0.0 24Aug11
* Version 3.0.0 Change to fix weights for simulated zero so mean overall is the same across the 2 groups
* Version 2.0.0 Seed is not mandatory
* Version 1.2.0 Fix to allow to nicely die if matrices are empty
* Version 1.1.0 
*		Also uses version 4 of QTE to deal with an issue with firpo
* 		reruns the firpo in sgkernel, then runs qte without firpo option	
*		Need to make sure pscore weights are recreated for each replicate
*** Below here is versions for sgnull.ado
* Version 1.0.0 
* 		Started from version 4.0.0 of sgnull
*               Allows subsampling to be with/without replacement
*               Removes own written bsample, uses canned
* 		Uses sample for without replacement	
* Version 4.0.0 adds subsampling option
*		Can select subsample size instead of full block BS
* Version 3.0.0 add zeros no time
*   		Fix something about bsreps/realdata realdata is incrementing by 1 each time, and bsreps = 1 for real data 0 all others
*
* Version 2.1.0: added time stuff, as well as fixing a few minor errors (which probly crashed the program)
*
* Version 1 never used, so no details there.
*

*** Taken from version 3 of sgnullsubsample
*** adds 2 outcomes which keep zeros in for figures with zeros

** CHANGELOG ** for sgnullsubsample
*!Version 3.0.0 14Dec11
*!Version 2.0.0 24Aug11
* Version 3.0.0 Change to fix weights for simulated zero so mean overall is the same across the 2 groups
* Version 2.0.0 Seed is not mandatory
* Version 1.2.0 Fix to allow to nicely die if matrices are empty
* Version 1.1.0 
*		Also uses version 4 of QTE to deal with an issue with firpo
* 		reruns the firpo in sgkernel, then runs qte without firpo option	
*		Need to make sure pscore weights are recreated for each replicate
*** Below here is versions for sgnull.ado
* Version 1.0.0 
* 		Started from version 4.0.0 of sgnull
*               Allows subsampling to be with/without replacement
*               Removes own written bsample, uses canned
* 		Uses sample for without replacement	
* Version 4.0.0 adds subsampling option
*		Can select subsample size instead of full block BS
* Version 3.0.0 add zeros no time
*   		Fix something about bsreps/realdata realdata is incrementing by 1 each time, and bsreps = 1 for real data 0 all others
*
* Version 2.1.0: added time stuff, as well as fixing a few minor errors (which probly crashed the program)
*
* Version 1 never used, so no details there.
*
*********************************************************************************************************


program define sgnullsubsampleznoweight, eclass sortpreserve

	syntax varlist (min=1 max=1) [if] [in] [aweight fweight /] , /*
		*/ IDvar(string)  Quantiles(string) SUBGroups(varlist) SAVing(string) TIMEvar(string) Treatmnt(varlist min=1 max=1)/*
		*/ [ BSreps(real 999) EVery(string) /*
		*/   SEED(string) model(string) covariates(varlist numeric) firpo(string asis) PSCOREname(varlist min=1 max=1)  Level(real 0.10) NSUBSAMP(string) REPLACEMENT(string)]

	preserve
	marksample touse
	qui keep if `touse'

	set matsize 2000
	if "`seed'" !="" {
		set seed `seed'
		}

	qui count
	local n = _result(1)

	tempvar weightvar
	if "`weight'" ~= "" { 
		gen double `weightvar'=`exp'
		local weight "[`weight'=`weightvar']"
	}
	else {
		gen double `weightvar'=1
		local weight "[aweight=`weightvar']"
	}
	

     	*parse the saving text
        parse `"`saving'"', parse(",")
        local saving "`1'"
        local savingcomma "`2'"
        local replace "`3'"

	*getting names of vars to be saved in results file into local macros savenames*
	local savenames_nt_nz "realdata"
	local savenames_nt_z "realdata"
	local savenames_nt_z_wz "realdata"
	local savenames_t_nz   "realdata"
	local savenames_t_z     "realdata"
	local savenames_t_z_wz     "realdata"
	forvalues q = `quantiles' { /* loop over quantiles numlist */
		/* note that use nz for nzs and z for zeros */
		/* 12/11 also note, saving an NZ version of t_z and nt_z for graphs */

		local savenames_nt_nz "`savenames_nt_nz' double yhq0nz_`q' double yhq1nz_`q' double ytq1_nt_nz_`q'"
		local savenames_nt_z   "`savenames_nt_z' double yhq0z_`q' double yhq1z_`q' double ytq1_nt_z_`q'"
		local savenames_nt_z_wz   "`savenames_nt_z_wz' double yhq0nz_`q' double yhq1nz_`q' double ytq1_nt_z_wz_`q'"
		local savenames_t_nz   "`savenames_t_nz' double yhq0nz_`q' double yhq1nz_`q' double ytq1_t_nz_`q'"
		local savenames_t_z     "`savenames_t_z' double yhq0z_`q' double yhq1z_`q' double ytq1_t_z`q'"
		local savenames_t_z_wz     "`savenames_t_z_wz' double yhq0nz_`q' double yhq1nz_`q' double ytq1_t_z_wz_`q'"

	} /* done with loop over quantiles numlist */

       	*set up post stuff by opening file handle
       	tempname bsresults_nt_nz bsresults_nt_z bsresults_nt_z_wz bsresults_t_nz bsresults_t_z bsresults_t_z_wz 
       	postfile `bsresults_nt_nz' bsrep `savenames_nt_nz' using `saving'_nt_nz `savingcomma' `replace' `every'
       	postfile `bsresults_nt_z'   bsrep `savenames_nt_z'   using `saving'_nt_z   `savingcomma' `replace' `every'
       	postfile `bsresults_nt_z_wz'   bsrep `savenames_nt_z_wz'   using `saving'_nt_z_wz   `savingcomma' `replace' `every'
       	postfile `bsresults_t_nz'   bsrep `savenames_t_nz'   using `saving'_t_nz   `savingcomma' `replace' `every'
       	postfile `bsresults_t_z'     bsrep `savenames_t_z'     using `saving'_t_z     `savingcomma' `replace' `every'
       	postfile `bsresults_t_z_wz'     bsrep `savenames_t_z_wz'     using `saving'_t_z_wz     `savingcomma' `replace' `every'


	*homemade preserve
	tempfile preserved_data
	qui save `preserved_data'
	forvalues bsrep=0/`bsreps' { /* BS loop */
		if `bsrep'>0 {
			di "." _continue
		}
		else {
			}

		local poststuff_nt_nz		/* clear last rep's results */
		local poststuff_nt_z		/* clear last rep's results */
		local poststuff_nt_z_wz		/* clear last rep's results */
		local poststuff_t_nz		/* clear last rep's results */
		local poststuff_t_z		/* clear last rep's results */
		local poststuff_t_z_wz		/* clear last rep's results */

		if `bsrep'>0 {			 /* bsample block */
			* for debug
			di c(seed)
			qui use `preserved_data', clear
			/* With replacement */
			if "`replacement'" ~="N" & "`replacement'" ~="n" {
				/* if using subsample */
				/* size of mini data set */
				if "`nsubsamp'" != "" {
					local N = `nsubsamp'
					*di "`N' is N"
					qui bsample `N', cluster(`idvar')
					}
				/* else use all data */
				else {
					qui bsample, cluster(`idvar')
					}
			} /* end with replacement block */
			/* Without replacement */
			else {
				/* Not subsample */
				if "`nsubsamp'" == "" {
					qui su `varlist' if `timevar'==1
					local N = _result(1)
					}
				/* Subsample */
				else {
					local N = `nsubsamp'
					}
				*di "N: `N'"
				qui sample `N' if `timevar'==1, count
				** keep only if first quarter there
				tempvar first_t idin
				qui gen `first_t' = `timevar'==1
			        qui egen `idin' = max(`first_t'), by(`idvar')
				qui keep if `idin'==1
		     	} /* end without replacement block */	
		} /* end bsample block */					

		if "`firpo'" != "" {
			* call sgkernel with firpo
			sgkernel `varlist' `weight', treatmnt(`treatmnt') subgroups(`subgroups') quantiles(`quantiles') timevar(`timevar') firpo(`firpo')

		}
		else    {
			sgkernel `varlist' `weight', treatmnt(`treatmnt') subgroups(`subgroups') quantiles(`quantiles') timevar(`timevar') 
		} /* end firpo block */

		tempname b 
		mat `b' = e(b)
		if `bsreps'==0 {
			tempname bpost v
			mat `bpost'=e(b)
			mat `v'=e(V)
			ereturn post `bpost' `v'

			di _new " -> Results for real data:" _new
			ereturn display
		}	


		forvalues q = `quantiles' { /* loop over quantiles numlist */

			tempname temphat0nz temphat0z temphat1nz temphat1z temptilde_nt_nz temptilde_nt_z temptilde_nt_z_wz temptilde_t_nz temptilde_t_z temptilde_t_z_wz

			mat `temphat0nz'  = `b'[1,"yhq0nz:q`q'"]	/* annoying stata rules require we save the value in a matrix*/
			local valuehat0nz = `temphat0nz'[1,1]		/* now we can save to a local macro */

			mat `temphat1nz'  = `b'[1,"yhq1nz:q`q'"]
			local valuehat1nz = `temphat1nz'[1,1]    

			mat `temphat0z'  = `b'[1,"yhq0z:q`q'"]
			local valuehat0z = `temphat0z'[1,1]	

			mat `temphat1z'  = `b'[1,"yhq1z:q`q'"]
			local valuehat1z = `temphat1z'[1,1]    

			mat `temptilde_nt_nz' = `b'[1,"ytq1_nt_nz:q`q'"]
			local valuetilde_nt_nz = `temptilde_nt_nz'[1,1]	

			mat `temptilde_nt_z' = `b'[1,"ytq1_nt_z:q`q'"]
			local valuetilde_nt_z = `temptilde_nt_z'[1,1]	

			mat `temptilde_nt_z_wz' = `b'[1,"ytq1_nt_z_wz:q`q'"]
			local valuetilde_nt_z_wz = `temptilde_nt_z_wz'[1,1]	

			mat `temptilde_t_nz' = `b'[1,"ytq1_t_nz:q`q'"]
			local valuetilde_t_nz = `temptilde_t_nz'[1,1]	

			mat `temptilde_t_z' = `b'[1,"ytq1_t_z:q`q'"]
			local valuetilde_t_z = `temptilde_t_z'[1,1]	

			mat `temptilde_t_z_wz' = `b'[1,"ytq1_t_z_wz:q`q'"]
			local valuetilde_t_z_wz = `temptilde_t_z_wz'[1,1]	

			mat drop `temphat0' `temphat1' `temptilde_nt_nz' `temptilde_nt_z' `temptilde_nt_z_wz' `temptilde_t_nz' `temptilde_t_z' `temptilde_t_z_wz'

			/* note here saving NZ versions with 0s include t_z_wz and nt_z_wz for graphs */
			local poststuff_nt_nz "`poststuff_nt_nz' (`valuehat0nz') (`valuehat1nz') (`valuetilde_nt_nz')"
			local poststuff_nt_z   "`poststuff_nt_z' (`valuehat0z') (`valuehat1z') (`valuetilde_nt_z')"
			local poststuff_nt_z_wz   "`poststuff_nt_z_wz' (`valuehat0nz') (`valuehat1nz') (`valuetilde_nt_z_wz')"
			local poststuff_t_nz   "`poststuff_t_nz' (`valuehat0nz') (`valuehat1nz') (`valuetilde_t_nz')"
			local poststuff_t_z     "`poststuff_t_z' (`valuehat0z') (`valuehat1z') (`valuetilde_t_z')"
			local poststuff_t_z_wz     "`poststuff_t_z_wz' (`valuehat0nz') (`valuehat1nz') (`valuetilde_t_z_wz')"
		} /* done with loop over quantiles numlist */							
	

		local realdata 0
	
		if `bsrep'==0 {
			local realdata = 1 
			di "setting realdata"
		}
		else {
		}

		post `bsresults_nt_nz' (`bsrep') (`realdata') `poststuff_nt_nz'
		post `bsresults_nt_z'   (`bsrep') (`realdata') `poststuff_nt_z'
		post `bsresults_nt_z_wz'   (`bsrep') (`realdata') `poststuff_nt_z_wz'
		post `bsresults_t_nz'   (`bsrep') (`realdata') `poststuff_t_nz'
		post `bsresults_t_z'     (`bsrep') (`realdata') `poststuff_t_z'
		post `bsresults_t_z_wz'     (`bsrep') (`realdata') `poststuff_t_z_wz'

	} /* done with BS loop */	

	postclose `bsresults_nt_nz'
	postclose `bsresults_nt_z'
	postclose `bsresults_nt_z_wz'
	postclose `bsresults_t_nz'
	postclose `bsresults_t_z'
	postclose `bsresults_t_z_wz'

	restore
end



program define sgkernel, sortpreserve

	syntax varlist(min=1 max=1) [if] [in] [aweight fweight /] , Quantiles(numlist) SUBGroups(varlist) TREATmnt(varlist min=1 max=1) TIMEvar(string) /*
		*/ [ model(string) covariates(varlist numeric) firpo(string asis) PSCOREname(varlist min=1 max=1) Level(string) ]

	preserve 
	marksample touse
	keep if `touse'


	if "`model'" ~= "" {
		local model "model(`model')"
       	}

	if "`covariates'" ~= "" {
		local covariates "covariates(`covariates')"
       	}

	if "`pscorename'" ~= "" {
		local pscorename "pscorename(`pscorename')"
       	}

	local depvar : word 1 of `varlist'

	tempvar weightvar
	/** change following to incorporate outcome of firpo if want */
	gen double `weightvar'=`exp'	/* we know that weight is nonmissing */


	/* run firpo if asked */
	if "`firpo'" ~= "" {
		di "in firpo" 
		tempvar pscore pscorewt
		`firpo' /* runs the model, which has to be fully specified by the user. compound quotes bc of asis option */
		qui predict double `pscore' if `touse'
		if "`pscorename'"~="" { 
			qui gen double `pscorename' = `pscore' if `touse' 
			}
		qui gen double `pscorewt' = `treatmnt'/`pscore' + (1-`treatmnt')/(1-`pscore') if `touse'
		su `pscorewt'
		/*  update weight */
		qui replace `weightvar'= `exp' * `pscorewt'
		su `weightvar'
	      	} /* end of firpo */


	/* Estimate mean tfx by subgroup */
	tempvar sg_var
	egen `sg_var' = group(`subgroups')	/* var that identifies subgroup membership uniquely */
	tab `sg_var' 
	local num_sg = _result(2)		/* number of subgroups */

	** define variables
	tempvar deltahat yt1_nt_nz yt1_nt_z yt1_nt_z_wz yt1_t_nz yt1_t_z yt1_t_z_wz nonzero newweight newweight2
	gen double `deltahat'       =.
	gen double `yt1_nt_nz'      = `depvar'    if `treatmnt'==1 	/* initialize */
	gen double `yt1_nt_z'       = `depvar'    if `treatmnt'==1 	/* initialize */

	*** NT_Z version with zeros in there
	gen double `yt1_nt_z_wz'    = `depvar'    if `treatmnt'==1 	/* initialize */
	gen double `yt1_t_nz'       = `depvar'    if `treatmnt'==1 	/* initialize */
	gen double `yt1_t_z'        = `depvar'    if `treatmnt'==1 	/* initialize */

	*** T_Z version with zeros in there
	gen double `yt1_t_z_wz'     = `depvar'    if `treatmnt'==1 	/* initialize */
	gen byte   `nonzero'    = `depvar'>0  

	*** weights to deal with zero adjustment
	gen double `newweight'  = `weightvar' if `treatmnt'==1
	gen double `newweight2' = `weightvar' if `treatmnt'==1

	*timevar stuff
	qui sum `timevar' if `touse'
	local tmin = _result(5)
	local tmax = _result(6)

	forvalues s=1/`num_sg' { /* loop over subgroups */

		*version without accounting for time or zeros
		* change to weightvar here not EXP so get right thing with firpo
		regress `depvar' `treatmnt' [`weight'=`weightvar'] if `sg_var'==`s' /* calculate mean treatment effect for this subgroup */
		replace `yt1_nt_nz' = `depvar' + _b[`treatmnt'] if `treatmnt'==0 & `sg_var'==`s' & e(sample)
				/* construct consis est of treated depvar for CG members under null */

		*version without accounting for time with zeros treated differently
		regress `depvar' `treatmnt' [`weight'=`weightvar'] if `sg_var'==`s' & `nonzero'==1 /* calculate mean tft for this subgroup given nonzero */
		replace `yt1_nt_z' = (`depvar' + _b[`treatmnt']) if `treatmnt'==0 & `sg_var'==`s' & `nonzero'==1
		replace `yt1_nt_z' =  0 		 	 if `treatmnt'==0 & `sg_var'==`s' & `nonzero'==0

		*Estimate nonzero shares for the null distribution, treating zeros DIFFERENTLY
		regress `nonzero' `treatmnt' [`weight'=`weightvar'] if `sg_var'==`s'
		local nonzero_share_`s'_0 = _b[_cons]			/* this is control group */
		local nonzero_share_`s'_1 = _b[_cons] + _b[`treatmnt']	/* this is tg            */

		*** first adjustment for nt_z, do same share of zeros with reweighting
		replace `newweight2'=`weightvar'*( /*
					*/   ( `nonzero_share_`s'_1' / `nonzero_share_`s'_0' ) * `nonzero' /*
					*/   + ( (1-`nonzero_share_`s'_1') / (1-`nonzero_share_`s'_0') ) * (1-`nonzero') /*
					*/ ) if `sg_var'==`s' & `treatmnt'==0 

		*varying by time 
		forvalues t=`tmin'/`tmax' { /* loop over time variable */
			*version with zeros treated like other obs
			regress `depvar' `treatmnt' [`weight'=`weightvar'] if `sg_var'==`s' & `timevar'==`t' /* calculate mean treatment effect for this subgroup */
			replace `yt1_t_nz' = `depvar' + _b[`treatmnt'] if `treatmnt'==0 & `sg_var'==`s' & `timevar'==`t' 
					/* construct consis est of treated depvar for CG members under null */

			*version with zeros treated differently
			regress `depvar' `treatmnt' [`weight'=`weightvar'] if `sg_var'==`s' & `timevar'==`t' & `nonzero'==1 /* calculate mean tft for this subgroup given nonzero */
			replace `yt1_t_z' = (`depvar' + _b[`treatmnt']) if `treatmnt'==0 & `sg_var'==`s' & `timevar'==`t' & `nonzero'==1
			replace `yt1_t_z' =  0 				if `treatmnt'==0 & `sg_var'==`s' & `timevar'==`t' & `nonzero'==0

			*Estimate nonzero shares for the null distribution, treating zeros DIFFERENTLY
			regress `nonzero' `treatmnt' [`weight'=`weightvar'] if `sg_var'==`s' & `timevar'==`t'
			local nonzero_share_`s'_`t'_0 = _b[_cons]			/* this is control group */
			local nonzero_share_`s'_`t'_1 = _b[_cons] + _b[`treatmnt']	/* this is tg            */


			* first adjustment for t_z, do same share of zeros with reweighting
			replace `newweight'=`weightvar'*( /*
							*/   ( `nonzero_share_`s'_`t'_1' / `nonzero_share_`s'_`t'_0' ) * `nonzero' /*
							*/   + ( (1-`nonzero_share_`s'_`t'_1') / (1-`nonzero_share_`s'_`t'_0') ) * (1-`nonzero') /*
							*/ ) if `sg_var'==`s' & `timevar'==`t' & `treatmnt'==0 

		} /* end loop over time variable */

	} /* done with loop over subgroups */
	

	* Initialize version of zero adjusted variables for later with zeros in 
	replace `yt1_nt_z_wz' = `yt1_nt_z'
	replace `yt1_t_z_wz' = `yt1_t_z'


	***********************************************************************************
	*** start, Chunk of code for debugging, checks new weights make			***
	*** sure that means are same within SG and overall for `depvar' and `nonzero'	***
	***********************************************************************************
	/* start of comment out debugging */
	* / * 
	* check first period/subgroup  means are the same

	di "Here first period/subgroup summarize nonzero shares"
	di "Here: nonzero `nonzero with weight weightvar `weightvar' real treatment sgvar 1"
	su `nonzero' if `treatmnt'==1 & `sg_var'==1 [aw=`weightvar']
	di "Here: nonzero `nonzero with weight newweight2 `newweight2' control sgvar 1"
	su `nonzero' if `treatmnt'==0 & `sg_var'==1 [aw=`newweight2']

	di "Here: nonzero `nonzero' with weight weightvar `weightvar' real treatment sgvar 1 time period 1"
	su `nonzero' if `treatmnt'==1 & `sg_var'==1 & `timevar' ==1 [aw=`weightvar']
	di "Here: nonzero `nonzero' with weight newweight `newweight' control time `timevar' sgvar 1 time period 1"
	su `nonzero' if `treatmnt'==0 & `sg_var'==1 & `timevar'==1 [aw=`newweight']

	di ""
	di ""
	di "Check means are same for real treatment and simulated ones in first sgvar subgroup"
	di ""
	di " ernq weightvar treatment"
	di "Command is: su `depvar' if `treatmnt'==1 & `sg_var' ==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 & `sg_var' ==1 [aw=`weightvar']
	di ""
	di "treatment `yt1_nt_nz' is yt1_nt_nz not conditional on positive, weight weightvar"
	di "Command is: su `yt1_nt_nz' if `treatmnt'==0 & `sg_var' ==1 [aw=`weightvar']"
	su `yt1_nt_nz' if `treatmnt'==0 & `sg_var' ==1 [aw=`weightvar']
	di ""

	di ""
	di ""
	di "Here: nt_z and ernq same means positive earnings (note that rest of values should be 0 so this should ensure same means)"
	di "check means are same for real treatment and yt1_nt_z first subgroup, positive"
	di "depvar `depvar' weight `weightvar' newweight newweight2"
	di ""
	di "command is: su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 [aw=`weightvar']
	di ""
	di "command is: su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 [aw=`newweight']"
	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 [aw=`newweight']
	di ""
	di "command is: su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1  [aw=`newweight2']"
	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1  [aw=`newweight2']
	di ""
	di "`yt1_nt_z' is yt1_nt_z given positive, weight `newweight2' newweight2"
	di "command is: su `yt1_nt_z' if `treatmnt'==0 & `depvar' > 0 & `sg_var' ==1 [aw=`newweight2']"
	su `yt1_nt_z' if `treatmnt'==0 & `depvar' > 0 & `sg_var' ==1 [aw=`newweight2']


	di ""
	di ""
	di "Here: nt_z and ernq same means overall "
	di "check means are ABOUT same for real treatment and yt1_nt_z first subgroup, overall"
	di "weightvar, newweight, newweight2 for treatment group"
	di ""
	di "command is: su `depvar' if `treatmnt'==1 &  `sg_var' ==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 &  `sg_var' ==1 [aw=`weightvar']
	di ""
	di "command is: su `depvar' if `treatmnt'==1 &  `sg_var' ==1 [aw=`newweight']"
	su `depvar' if `treatmnt'==1 &  `sg_var' ==1 [aw=`newweight']
	di ""
	di "command is: su `depvar' if `treatmnt'==1 &  `sg_var' ==1  [aw=`newweight2']"
	su `depvar' if `treatmnt'==1 &  `sg_var' ==1  [aw=`newweight2']
	di ""
	di "`yt1_nt_z' is yt1_nt_z overall, weight `newweight2' newweight2"
	di "command is: su `yt1_nt_z' if `treatmnt'==0 &  `sg_var' ==1 [aw=`newweight2']"
	su `yt1_nt_z' if `treatmnt'==0 &  `sg_var' ==1 [aw=`newweight2']

	di ""
	di ""
	di "t_nz and ernq same means overall"
	di ""
	di "ernq first time period and sgvar weightvar"
	di "command is: su `depvar' if `treatmnt'==1 & `sg_var' ==1 & `timevar'==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 & `sg_var' ==1 & `timevar'==1 [aw=`weightvar']
	di ""
	di "`yt1_t_nz' is yt1_t_nz overall, weight `weightvar' weightvar"
	di "command is: su `yt1_t_nz' if `treatmnt'==0 & `sg_var' ==1 & `timevar'==1 [aw=`weightvar']"
	su `yt1_t_nz' if `treatmnt'==0 & `sg_var' ==1 & `timevar'==1 [aw=`weightvar']

	di ""
	di "t_z and ernq same means positive"
	di ""
	di "command is:	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`weightvar'] 
	di ""
	di "command is:	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`newweight']"
 	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`newweight']
	di ""
	di "command is:	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`newweight2']"
	su `depvar' if `treatmnt'==1 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`newweight2']
	di ""
	di "`yt1_t_z' is yt1_t_z conditional on positive, , weight `newweight' newweight"
	di "command is: su `yt1_t_z' if `treatmnt'==0 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`newweight']"
	su `yt1_t_z' if `treatmnt'==0 & `depvar' > 0 & `sg_var' ==1 & `timevar'==1 [aw=`newweight']


	di ""
	di "t_z and ernq same means period 1"
	di ""
	di "command is:	su `depvar' if `treatmnt'==1  & `sg_var' ==1 & `timevar'==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1  & `sg_var' ==1 & `timevar'==1 [aw=`weightvar'] 
	di ""
	di "command is:	su `depvar' if `treatmnt'==1  & `sg_var' ==1 & `timevar'==1 [aw=`newweight']"
 	su `depvar' if `treatmnt'==1  & `sg_var' ==1 & `timevar'==1 [aw=`newweight']
	di ""
	di "command is:	su `depvar' if `treatmnt'==1  & `sg_var' ==1 & `timevar'==1 [aw=`newweight2']"
	su `depvar' if `treatmnt'==1  & `sg_var' ==1 & `timevar'==1 [aw=`newweight2']
	di ""
	di "`yt1_t_z' is yt1_t_z overall, , weight `newweight' weight"
	di "command is: su `yt1_t_z' if `treatmnt'==0  & `sg_var' ==1 & `timevar'==1 [aw=`newweight']"
	su `yt1_t_z' if `treatmnt'==0  & `sg_var' ==1 & `timevar'==1 [aw=`newweight']
	pause 

	* check overall means are about the same
	* noi di "HERE overall means"

	di "Here: `nonzero with weight `weightvar' real treatment full sample"
	di ""
	di "Command is:	su `nonzero' if `treatmnt'==1 [aw=`weightvar']"
	su `nonzero' if `treatmnt'==1 [aw=`weightvar']
	di ""
	di "`nonzero with weight `newweight' control full sample"
	di "Command is:	su `nonzero' if `treatmnt'==0 [aw=`newweight']"
	su `nonzero' if `treatmnt'==0 [aw=`newweight']
	di ""
	di "`nonzero with weight `newweight2' control full sample"
	di "Command is:	su `nonzero' if `treatmnt'==0 [aw=`newweight2']"
	su `nonzero' if `treatmnt'==0 [aw=`newweight2']
	di ""
	di ""
	di "Check overall means about equal"
	di ""
	di ""
	di ""
	di "Earnings in Treatment group, weightvar"
	di "Command is:	su `depvar' if `treatmnt'==1 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 [aw=`weightvar']
	di ""
	di "yt1_nt_nz, weightvar, control"
	di "Command is: su `yt1_nt_nz' if `treatmnt'==0 [aw=`weightvar']"
	su `yt1_nt_nz' if `treatmnt'==0 [aw=`weightvar']
	di ""
	di "yt1_t_nz, weightvar, control"
	di "Command is: su `yt1_t_nz' if `treatmnt'==0 [aw=`weightvar']"
	su `yt1_t_nz' if `treatmnt'==0 [aw=`weightvar']
	di ""
	di "yt1_nt_z, newweight2, control"
	di "Command is: su `yt1_nt_z' if `treatmnt'==0 [aw=`newweight2']"
	su `yt1_nt_z' if `treatmnt'==0 [aw=`newweight2']
	di ""
	di "yt1_nt_z, newweight2, control"
	di "Command is: su `yt1_nt_z' if `treatmnt'==0 [aw=`newweight2']"
	su `yt1_nt_z' if `treatmnt'==0 [aw=`newweight2']

	di ""
	di ""
	di " Check means positive zeros"
	di ""
	di "Earnings in treatment group, weightvar, positive"
	di "Command is:	su `depvar' if `treatmnt'==1 & `depvar' > 0 [aw=`weightvar']"
	su `depvar' if `treatmnt'==1 & `depvar' > 0 [aw=`weightvar']
	di ""
	di "yt1_nt_z newweight2, control, positive"
	di "Command is:	su `yt1_nt_z' if `treatmnt'==0 & `depvar' > 0 [aw=`newweight2']"
	su `yt1_nt_z' if `treatmnt'==0 & `depvar' > 0 [aw=`newweight2']
	di ""
	di "yt1_t_z newweight, control, positive"
	di "Command is:	su `yt1_t_z' if `treatmnt'==0 & `depvar' > 0 [aw=`newweight']"
	su `yt1_t_z' if `treatmnt'==0 & `depvar' > 0 [aw=`newweight']
	pause	

	***********************************************************************************
	*** end, Chunk of code for debugging, checks new weights make			***
	*** sure that means are same within SG and overall for `depvar' and `nonzero'	***
	***********************************************************************************
	*  / /* end of comment out debugging */

	
	* This Version has different yhq0 and yhq1 when treating zeros differently
	* yhq0nz and yhq1nz are old version, treating zeros the same
	* yhq0z and yhq1z are zero different version, estimated only on nonzero obs
	tempname b yhq0nz yhq0z yhq1nz yhq1z ytq1_nt_nz ytq1_nt_z ytq1_nt_z_wz ytq1_t_nz ytq1_t_z ytq1_t_z_wz beta

	*Estimate real-data qte TREATING ZEROS THE SAME 
	*** now using weightvar, no firpo (firpo is used to create weightvar above)
	di "1 Real Data Zeros in: qte `depvar' `treatmnt' [`weight'=`weightvar'] , quantiles(`quantiles') nokd nojcumul "
       	qte `depvar' `treatmnt' [`weight'=`weightvar'] , quantiles(`quantiles') nokd nojcumul 	/* hardcoded options: nokd, nojcumul no firpo, weights with firpo above*/
	mat `b' = e(b)

	mat `yhq0nz' = `b'[1,"c:"]	/* this is vector of quantiles for actual cg */
	mat coleq `yhq0nz' = yhq0nz	/* set names to yhq0nz for real-data tg quantiles */

	mat `yhq1nz' = `b'[1,"t:"]	/* this is vector of quantiles for actual tg */
	mat coleq `yhq1nz' = yhq1nz	/* set names to yhq1nz for real-data tg quantiles */

	*Estimate real-data qte TREATING ZEROS DIFFERENTLY, only do for nonzero obs
	di "2 Real Data Zeros Out: qte `depvar' `treatmnt' [`weight'=`weightvar']  if `nonzero' != 0, quantiles(`quantiles') nokd nojcumul "
       	qte `depvar' `treatmnt' [`weight'=`weightvar'] if `nonzero' != 0, quantiles(`quantiles') nokd nojcumul 	/* hardcoded options: nokd, nojcumul no firpo*/
	mat `b' = e(b)

	mat `yhq0z' = `b'[1,"c:"]	/* this is vector of quantiles for actual cg */
	mat coleq `yhq0z' = yhq0z	/* set names to yhq0z for real-data tg quantiles */

	mat `yhq1z' = `b'[1,"t:"]	/* this is vector of quantiles for actual tg */
	mat coleq `yhq1z' = yhq1z	/* set names to yhq1z for real-data tg quantiles */

	*Estimate qte for the null distribution, treating zeros LIKE other obs, nt_nz
	di "3 NT_NZ: qte `yt1_nt_nz' `treatmnt' [`weight'=`weightvar'] , quantiles(`quantiles') `model' `covariates'   `pscorename' nokd nojcumul "
	qte `yt1_nt_nz' `treatmnt' [`weight'=`weightvar'] , quantiles(`quantiles') `model' `covariates'   `pscorename' nokd nojcumul 
	/*notg*/ 	/* hardcoded options: nokd, nojcumul no firpo */

	mat `b' = e(b)
	mat `ytq1_nt_nz' = `b'[1,"c:"]			/* this is vector of synthetic quantiles for TG based on (actual CG plus mean treatment effect) */
	mat coleq `ytq1_nt_nz' = ytq1_nt_nz		/* set names to ytq for null tg quantiles */

	*Estimate qte for the null distribution, treating zeros DIFFERNTLY, nt_z, newweight2
	di "4 NT_Z, Zeros Out: qte `yt1_nt_z' `treatmnt' [aweight=`newweight2'] if `nonzero' != 0, q(`quantiles') nokd nojcumul notg"
	qte `yt1_nt_z' `treatmnt' [aweight=`newweight2'] if `nonzero' != 0, q(`quantiles') nokd nojcumul notg
	mat `b' = e(b)	/* already created tempname for b above and can now recycle it */

	mat `ytq1_nt_z' = `b'[1,"c:"]		/* this is vector of quantiles for TG based on (actual CG plus mean treatment effect) */
	mat coleq `ytq1_nt_z' = ytq1_nt_z	/* set names to ytq for null tg quantiles */

	*Estimate qte for the null distribution, treating LIKE other obs, nt_z_wz, newweight2
	di "4alt NT_Z, Zeros In:	qte `yt1_nt_z_wz' `treatmnt' [aweight=`newweight2'], q(`quantiles') nokd nojcumul notg"
	qte `yt1_nt_z_wz' `treatmnt' [aweight=`newweight2'], q(`quantiles') nokd nojcumul notg
	mat `b' = e(b)	/* already created tempname for b above and can now recycle it */

	mat `ytq1_nt_z_wz' = `b'[1,"c:"]		/* this is vector of quantiles for TG based on (actual CG plus mean treatment effect) */
	mat coleq `ytq1_nt_z_wz' = ytq1_nt_z_wz	/* set names to ytq for null tg quantiles */

	*Estimate qte for the null distribution, treating zeros LIKE other obs, t_nz
	di "5 T_NZ: qte `yt1_t_nz' `treatmnt' [`weight'=`weightvar'] , quantiles(`quantiles') nokd nojcumul notg  "
	qte `yt1_t_nz' `treatmnt' [`weight'=`weightvar'] , quantiles(`quantiles') nokd nojcumul notg	/* hardcoded options: nokd, nojcumul no firpo */
	mat `b' = e(b)	/* already created tempname for b above and can now recycle it */

	mat `ytq1_t_nz' = `b'[1,"c:"]			/* this is vector of synthetic quantiles for TG based on (actual CG plus mean treatment effect) */
	mat coleq `ytq1_t_nz' = ytq1_t_nz		/* set names to ytq for null tg quantiles */

	*Estimate qte for the null distribution, treating zeros DIFFERNTLY, t_z, newweight
	di "6 T_Z, Zeros Out: qte `yt1_t_z' `treatmnt' [aweight=`newweight']  if `nonzero' != 0, q(`quantiles') nokd nojcumul notg"
	qte `yt1_t_z' `treatmnt' [aweight=`newweight']  if `nonzero' != 0, q(`quantiles') nokd nojcumul notg

	mat `b' = e(b)	/* already created tempname for b above and can now recycle it */

	mat `ytq1_t_z' = `b'[1,"c:"]		/* this is vector of quantiles for TG based on (actual CG plus mean treatment effect) */
	mat coleq `ytq1_t_z' = ytq1_t_z		/* set names to ytq for null tg quantiles */

	*Estimate qte for the null distribution, treating zeros DIFFERNTLY, t_z_wz, newweight
	di "6alt, T_Z, Zeros In: qte `yt1_t_z_wz' `treatmnt' [aweight=`newweight'] , q(`quantiles') nokd nojcumul notg"
	qte `yt1_t_z_wz' `treatmnt' [aweight=`newweight'], q(`quantiles') nokd nojcumul notg

	mat `b' = e(b)	/* already created tempname for b above and can now recycle it */

	mat `ytq1_t_z_wz' = `b'[1,"c:"]		/* this is vector of quantiles for TG based on (actual CG plus mean treatment effect) */
	mat coleq `ytq1_t_z_wz' = ytq1_t_z_wz		/* set names to ytq for null tg quantiles */

	*create matrices for posting
	mat `b' = `yhq0nz' , `yhq1nz' , `yhq0z' , `yhq1z' , `ytq1_nt_nz' , `ytq1_nt_z' , `ytq1_nt_z_wz' , `ytq1_t_nz' , `ytq1_t_z' , `ytq1_t_z_wz' 

	*create dummy "variance" matrix for posting purposes
	tempname v

	mat `v' = J(colsof(`b'),colsof(`b'),0)

	local colnamesb: colfullnames(`b')
	mat colnames `v' = `colnamesb'
	mat rownames `v' = `colnamesb'
	
	* check that matrices are not missing
        tempvar y0nz y1nz y0z y1z ytntnz ytntz yttnz yttz ytntzwz yttzwz
	gen `y0nz' = matmissing(`yhq0nz')
	gen `y1nz' = matmissing(`yhq1nz')
	gen `y0z' = matmissing(`yhq0z')
	gen `y1z' = matmissing(`yhq1z')
	gen `ytntnz' = matmissing(`ytq1_nt_nz')
	gen `yttnz' = matmissing(`ytq1_t_nz')
	gen `ytntz' = matmissing(`ytq1_nt_z')
	gen `yttz' = matmissing(`ytq1_t_z')
	gen `ytntzwz' = matmissing(`ytq1_nt_z_wz')
	gen `yttzwz' = matmissing(`ytq1_t_z_wz')

	* if any are missing, save data, don't return results
	tempname crashed
	if `y0nz'==1 | `y1nz'==1 | `y0z'==1 | `y1z'==1 | `ytntnz'==1 | `ytntz'==1 | `yttnz'==1 | `yttz'==1 | `ytntzwz'==1 | `yttzwz'==1 {
		* check if crashed already
		di "`crashed'"
		capture confirm file "`saving'crashed"
		if _rc==0 {
			di c(seed)
			append using `saving'crashed
			save `saving'crashed, replace
			}
		else 	{
			save `saving'crashed, replace
			exit
			}
		}
	else {
		di "firpo `firpo'"
		*post results
		ereturn post `b' `v' , depname(`depvar') esample(`touse')
		*ereturn display
		}
end

